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ABSTRACT 

We consider preconditioning methods to accelerate convergence to a steady state for the 
incompressible fluid dynamic equations. The analysis relies on the inviscid equations. The 
preconditioning consists of a matrix multiplying the time derivatives. Thus the steady state 
of the preconditioned system is the same as the steady state of the original system. We 
compare our method to other types of pseudo-compressibility. For finite difference methods 
preconditioning can change and improve the steady state solutions. An application to viscous 
flow around a cascade with a non-periodic mesh is presented. 
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1 Introduction 


One way to solve the steady state incompressible equations is to march the time dependent equations 
until a steady state is reached. Since the transient is not of any interest one can use acceleration 
techniques which destroy the time accuracy but enables one to reach the steady state faster. Such 
methods can be considered as preconditionings to accelerate the convergence to a steady state. 
For the incompresible equations the continuity equation does not contain any time derivatives. 
To overcome this difficulty, Chorin [3] added an artificial time derivative of the pressure to the 
continuity equation together with a multiplicative variable, (3 . With this artificial term the resultant 
scheme is a symmetric hyperbolic system for the inviscid terms. Thus, the system is well posed 
and and numerical method for hyperbolic systems can be used to advance this system in time.. 
The free parameter (3 is then chosen to reach the steady state quickly. Later Turkel ([8], [9], [10]) 
extended this concept by adding a pressure time derivative also to the momentum equations. The 
resulting system after preconditioning is no longer symmetric but can be symmetrized by a change 
of variables. 

Thus, we will consider systems of the form 

4* fx 4* 9y — 0 * 

This system is written in conservation form though for some applications this is not necessary. Our 
analysis will be based on the linearized equations so the conservation form does not appear in the 
analysis though it does appear in the final numerical approximation. This system is now replaced 

hy 


P 1 v) t + f x 4- g y = 0, (1) 

or in linearized form 

P -1 w ( + Aw x 4- Bw y = 0, (2) 

with A and B constant matrices. 

For this system to be equivalent to the original system, in the steady state, we demand that 
P -1 have an inverse. This only need be true in the flow regime under consideration. We shall 
see later that frequently P is singular at stagnation points. Thus, we will temporarily consider 
strictly flows without a stagnation point. We also assume that the Jacobian matrices A = |£ 
and B = are simultaneously symmetrizable. In terms of the ‘symmetrizing’ variables we also 
demand that P be positive definite. We shall show later, in detail, that it does not matter which 
set of dependent variables are used to develop the preconditioner. One can transform between any 
two sets of variables. Thus, when we are finished we will analyze a system which is similar to (2), 
where the matrices A and B are symmetric and P is both symmetric and positive definite. Such 
systems are known as symmetric hyperbolic systems. One can then multiply this system by w and 
integrate by parts to get estimates for the integral of wf, i.e. energy estimates. These estimates 
can then be used to show that the system is well posed . We stress that if P is not positive then 
we may change the physics of the problem. For example, if P = —I then we have reversed the time 
direction and must therefore change all the boundary conditions. Keeping the right signs for the 
eigenvalues is a necessary but not sufficient condition for well-posedness. 

With this assumption the steady state solutions of the two systems are the same. Assuming 
that the steady state has a unique solution, it does not matter which system we march to a steady 
state. We shall later see that for the finite difference approximations the steady state solutions are 
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not necessarily the same and usually the preconditioned system leads to a better behaved steady 
state. 

2 Incompressible Equations 

Consider the incompressible inviscid equations in primitive variables. 


Wj ~l” Vy 0 

U t + UU X + VUy +/>* = () 

Vi + UV X + VVy + Py = 0 

We generalize Chorin’s pseudo-compressibility method [3]. Using the preconditioning suggested in 
[8] (with a — 1) we have 


1 

02 Pt + U X + Vy = 0 

u 

—p t + U t + uu x + VUy + p x = 0 

V 

■jpPt + V t + UV X -\ - VVy + Py - 0 (3) 

or in conservation form 

1 

02 Pt + Ux + V y 

2 'll 

-0^Pt + U t + (u 2 + p) x + (uv)y 
2v 

02Pt + V t + (uv) x + (v 2 + p) y 
We can also write (3) in matrix form using 

/ 1//3 2 0 0 \ 

p - 1 = uta 2 io, 

\ v/P 2 0 1 j 

/ P 2 0 0 \ 

P = -« 1 0 

\ -v 0 l ) 

Multiplying by P we rewrite this as 

w t + PAw x + PBwy = 0. 

Let 

D = u>i A + u> 2 B - 1 < ui,W2 < 1 


= 0 
= 0 
= 0 
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where are the Fourier transform variables in the x and y directions respectively. The speeds 

of the waves are now governed by the roots of det(\I — P Au>\ — P Buz) = 0 or equivalently 
det(XP * — Auj\ — Buj 2 ) = 0. Let 


q — uu 1 + vu>2 . 


Then the eigenvalues of PD are 


do = q ( 4 ) 

d± = ±(3 

and so the ‘acoustic’ speed is isotropic. 

The spatial derivatives involve symmetric matrices, i.e. D is a symmetric matrix but P is not 
symmetric. Thus, while the original system was symmetric hyperbolic the preconditioned system 
is no longer symmetric. In [8] it is shown that as long as 

fl 2 > (u 2 4- v 2 ) 

the equations can be symmetrized. On the other hand the eigenvalues are most equalized if /3 2 = 
(u 2 + v 2 ) [8]. So, we wish to choose (3 2 slightly larger than u 2 + v 2 . However, numerous calculations 
verify, that in general, a constant (3 is the best for the convergence rate. The reasons for this are 
not clear. 

We wish to stress that (3 has the dimensions of a speed. Therefore, /? cannot be a universal 
constant. There are papers that claim that (3 = 1 or (3 = 2.5 are optimal. Such claims cannot be 
true in general. It is simple to see that if one nondimensionalizes the equation then (3 gets divided 
by a reference velocity. Hence, the optimal ‘constant’ /? depends on the dimensionalization of the 
problem and in particular depends on the inflow conditions. In many calculations the inflow mass 
flux is equal to 1 or alternatively p + (u 2 + v 2 )/2 = 1. Such conditions will give an optimal (3 close 
to one. 

We next define the Bernoulli function 

H =p + (u 2 + v 2 )/ 2. 

Bernoulli’s theorem states that when the flow is steady and inviscid then H is constant along 
streamlines. We now multiply the second equation of (3) by u and the third equation of (3) by v 
and add these two equations. If /3 2 = u 2 + v 2 , the result is 


Ht + uH x + vHy — 0. (5) 

Thus, by altering the time dependence of the equations we have constructed a new equation in 
which H is convected along streamlines. Furthermore, if H is a uniform constant both initially and 
at inflow then H will remain constant for all time. On the numerical level this will usually not be 
true because of the introduction of an artificial viscosity or because of upwinding. For viscous flow, 
(5) is replaced by 


H t + uH x + vH y = — (uAu + »A«) 
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We note that these relationships for H follow from the momentum equations and do not depend 
on the form of the continuity equation. Hence, we consider the following generalization of (3) 

jpPt + aH t + u x + Vj , = 0 


au 


-jpPt + + uu x -|- vu y + p x = 0 


av 


-jpPt + V t + UV X +VVy +Py = 0 

where, a, a and 0 are free parameters. When u> 2 + u;| = 1 the eigenvalues of PD are 

s ± v/« 2 + Wd 


( 6 ) 


UJ\U + U>2V, 


2d 


where 


Q 2 

q 2 = u 2 + v 2 , d = 1 + a — a^, s = (1 — a)((*qu + u^v). 


Hence, the ‘acoustic’ eigenvalue is isotropic if a = 1. Furthermore, d = 1 if either a = 0 or 
P 2 = u 2 + v 2 . For a = 0 we recover our original scheme. For a = -1 the time derivative of the 
pressure no longer appears in the continuity equation. For general a,0 we have 


P _1 = 4r 


a 2 


p = 


1 + a - aa^p- 


( (a 


au 

av \ 



au 

p 2 

0 



av 

0 

p 2 ) 


0 2 


—au 


—av 

—au 1 + 

a- 

w 

aauv 

p* 


\ 


-av 


aauv 


1 + a - ctav 


2 

w- ) 


If we write the equation in conservation form (1) we have 


>-i 

conservative 


1 

w 


( ( a +l) au av 

(1 + a + a)u fi 2 + au 2 auv 

(1 -f a + a)v auv /3 2 + av 2 


conservative — 


1 + a - aa^p- 


( P 2 + a(u 2 + v 2 ) — au 

— (1 + a + a)u i + a-*p 
i -(l + a + a)v 


—av 

a auv 


aauv 


1 + a - ^3L 


w ) 


In [9] an analogy to the symmetric preconditioning of van Leer, Lee and Roe was constructed 
for the incompressible equations. If we choose a = 1 P is symmetric. If we also choose 0 2 = u 2 + v 2 
then we get the preconditioning of van Leer et.al. . 


P = 


u 2 - \-v 2 — u 

l+ds 


—u 
- v 


—v 

uv 

uv i , y 2 

t t ~T~l 




u 2 +v 2 


u 2 +u 2 
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These examples show that preconditioning is not unique. If fact, since the determinant of the 
transpose of a matrix is equal to the determinant of the original matrix it follows that the transpose 
of P is also a preconditioner with the same eigenvalues for the preconditioned system. These various 
systems will have the same eigenvalues but different eigenvectors for the preconditioned system. 
Numerous calculations show that the system given by P in (3) is more robust and converges faster 
than with the transpose preconditioner. This shows that it is not sufficient to consider just the 
eigenvalues but that the eigenvectors are also of importance. The eigenvectors are given in ([10]). 
We next consider the preconditioner considered by de Jouette, Viviand et. al. ([4]). Define: 

q = uu x + vu y + p x - ( t xx + T xy ) 

r = UV X + Wy +Py - {T X y + Tyy) 

s = uq + vr 

U 2 = u 2 + v 2 

Then they consider the following extension of the incompressible Navier-Stokes equations. 


Pt + dy(u x + v y ) + ays = 0 

u t + ayq + + v y) + eyus = 0 (7) 

v t + ayr + (3yv(u x + v y ) + eyvs = 0 

In the steady state q = r = s = 0 and u x + v y = 0 and so we recover the usual incompressible 

equations, ay ,dy ,ey ,ay ,j3y are free parameters that satisfy the following conditions 


ay/3y — dy€y 
ay > 0 dy > 0 
(dy + a V U 2 )(dy + PyU 2 ) >0 

In addition, in order for the speed of the convective wave to remain unchanged we add the condition 
ay = 1. From the momentum equations we obtain 

[uu t + vv t + flyU 2 (u x + v y )] 

5__ 1 + eyU 2 


Hence we can rewrite (8) as 


1 + eyU 2 a v . . 

d Pt ~ dv + VVt ' 

+ U X +Vy= 0 

3yu „ 

- —j~Pt + u t + q = 0 

fivv 

- "5 Pt + v t + r = 0 

Comparing this with (6) we see that the two approaches are identical if 
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dy = 

0 * 

(a + l )/? 2 — aaU 2 

ay — 


Uv = 

~p iv 

ey = 

aa 

p dv 


Choosing a = 1 and /? = U 2 we get the standard preconditioning (3). The Viviand parameters 
become ay = —a, j3y = — 1, ay = 1, dy = U 2 , ey = a . Then a — 0 gives the Turkel preconditioner 
and a = 1 gives the van-Leer (symmetric) preconditioner. 


3 Difference equations 


Until now the entire analysis has been based on the partial differential equation. We now make 
some remarks on important points for any numerical approximation of this system. When using 
a scheme based on a Riemann solver this solver should be for the preconditioned system and not 
the original scheme. When using a central difference schemes there is a need to add an artificial 
viscosity. Accuracy is improved for low Mach number flows if the preconditioner is applied only 
to the physical convective and viscous terms but not to the artificial viscosity. The use of a 
matrix artificial dissipation ([7]) should be based on the preconditioned equations as for Riemann 
solvers, difference scheme. Hence, both for upwind and central difference schemes the Riemann 
solver or artificial viscosity should be based on P -1 |PA| and not |A| i.e. in one dimension solve 
Wt + P/x = P(P -1 |PA|tn r ) x . When using characteristics for extrapolation at the boundaries 
it should be based on the characteristics of the modified system and not the physical system. 
Preconditioning is even more important when using multigrid than with an explicit scheme. With 
the original system, the stiffness of the eigenvalues greatly affects the smoothing rates of the slow 
components and so slows down the multigrid method, [6]. We conclude that the steady state 
solution of the preconditioned system may be different from that of the physical system. Thus, on 
the finite difference level the preconditioning can improve the accuracy as well as the convergence 
rate. 

We next consider adding artificial viscosity to the system (3). We first rewrite this system 
eliminating p t from the velocity equations. This gives 


Pt + P 2 (u x + v y ) = 

U t +P X + VUy - UVy = 
- VU X + Py = 


or in matrix form 



v t + uv. 


( o /? 2 o \ 

+ l o o 

y o —v u J 


( p \ / 0 0 

u I + I 0 v 

W* V° 0 


0 

0 

0 



= 0 ( 8 ) 
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w t + P Aw x + P Bw y = 0. 


We next consider the use of a matrix valued viscosity. Let D = oj\A + u 2 B with u> 2 + u> 2 = 1. 
The non-preconditioned matrix viscosity is given by |j 4| in the x direction and |B| in the y direction. 

Then 


\D\ = 


1 

(a 2 -a 3 )<? 2 


2 q 2 
qR 2 
v qRS 


q R 2 

(a 2 2 +a 3 2 ) r 2 + (a 2 - A 3 ) S 2 1#| 

RS (A 2 2 + A 3 2 -(A 2 -A 3 )|fl|) 


qRS \ 

R S (a 2 2 + A 3 2 — (A 2 — A 3 ) \R\j 

(a 2 2 + a 3 2 ) S 2 + (A 2 - A 3 ) R 2 |fl| J 


R + y/W+4 , _ R-y/W+J 

A 2 - 2 , A 3 - 2 

R = UU) 1 + VU) 2 , S - UU 2 — 11 W 2 

For the preconditioned artificial viscosity we consider instead P -1 |PA| and P -1 |PB| (see [7]). 
We consider the case a — 1 with P and a arbitrary. Then 

/ V n V\2 V X3 \ 

P" 1 |PD| = V n V 22 V 23 . 

y V 31 Vi2 V 33 J 


with 

tr 1 + a aS 2 q 2 X 

v " = J7d + 


Via = a 


pVd 


+ S 2 uX + 


RS (f3 2 — q 2 ) vX 


v * = J7= d + s2uX + 


RS ( (3 2 -aq 2 ) vX 

ft 2 


V 13 = a 


7 = + SvX 35 

/?v/S /? 2 


„ a (/? 2 -a<? 2 ) 11 A 

V 31 = + 5 v X — a2 


P y/d 


p p 2 S 2 u 2 X 


- RS 1 + a - -V uvX + 


2 P 2 \ v fi 2 {jP-dq*} ti 2 X 


) 


33 sTd q 2 


+ RS 1 + 0 - + 


2 /3 2 \ v , R 2 (P 2 -d^)v 2 X 


) 
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where 


C 0 2 Sv+Ru (q 2 ~P 2 )) (f3 2 Su+ Rv(/3 2 -aq 2 )) X 

P 2 q 2 


T/ _(P 2 Sv + Ru (aq 2 — f3 2 )) (f3 2 Su + Rv {(5 2 - q 2 )) X 
32 /? 2 q 2 


UUq + VU> 2 UU>2 — VOJ i 

K = o = 


d = 1 + a - a^— 9 2 = « 2 + w 2 


X = 




dZ 


„ _ P 4 S 2 - (f3 2 - q 2 ) (f3 2 - aq 2 ) R 2 

z F 

By inspection the matrix is symmetric when a = 1 . For the special case a = 1 and (3 2 = u 2 + v 2 
the formulas simplify and we get 


v„ = ^ 

9 


V 21 = V 12 = R- 
9 

V 31 = Vj3 = R- 
9 

U 2 (ft-1) t> 2 + f?tt 2 

V22 = 9 + — = 

9 9 

v 2 (#-l) u 2 + rw 2 

^33 = 9 + — = — 

9 9 


V23 = V32 


uv(R — 1) 
9 


For the equations in conservation form we multiply the continuity equation by u and add to 
the x velocity equation. We also multiply the continuity equation by v and add to the y velocity 
equation. 
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4 Computational Results 

We now present a calculation for two dimensional flow around an cascade to demonstrate the 
previous theory. The discretization is based on the multistage time method coupled with a central 
difference approximation as described in ([5], [7]). The basic scheme is accelerated by using a local 
time step, residual smoothing and multigrid. This code was further developed to consider cascade 
configurations in which the grid is not necessarily continous across the wake ([1], [2])- We compute 
the flow about a NACA0012 with periodic external boundaries. The flow is turbulent and we use 
a Baldwin-Lomax turbulence model, with Re = 500,000, Pr = 0.7, Pr t = 0.9 At inflow the angle 
of attack is specified as well as the Bernoulli constant, p + - - 1. The mesh is 192 x 32 and 

is shown in figure 1. We use a four stage Runge-Kutta method as a smoother for a full multigrid 
iteration. We choose a = 0 and 0 2 3 = max(K(u 2 + v 2 ),0 m ) with A = 1.1 ,0 m = 0.4, see (6). In 
figure 2 we plot the convergence rate for different values of o. We see that the fastest convergence 
occurs when a = 1 followed by a = 0 and finally a ~ -1. We also considered viscous flow about 
a VKI cascade (figure 3). In this case the convergence of all the methods slowed down, a = 1 was 
still the most efficient method but the differences were less dramatic than in the previous case. In 
other cases in was necessary to choose 0 almost constant. The symmetric preconditioner, a 1 
was more robust but not faster than a = 0. 


5 Conclusions 

A three parameter preconditioning matrix has been introduced for the incompressible inviscid 
equations. This is equivalent to the pseudo- compressibility methods considered by de Jouette et. 
al. When a = 1 the ‘acoustic’ speeds are symmetric. Furthermore, one can choose the parameter 
a so that the preconditioning matrix is symmetric. For the inviscid case considered computed a 
considerable increase in the convergence rate was achieved. 

In addition the incompressible equations offer a theoretical advantage over the compressible 
equations for the theoretical study of preconditioning methods. This is because of the simpler 
nature of the equations and the fact that the original method of Chorin is already symmetric. 
Nevertheless, a central difference scheme coupled with a Runge-Kutta time advancement suffers 
from lack of robustness. In particular 0 needs to be bounded away from zero at a relatively 
high level for many of the cases. Using the symmetric preconditioner a = a = 1 yields a more 
robust scheme though it does not seem to converge faster than the nonsymmetric preconditioner. 
Furthermore, changes of the physical inflow boundary condition can greatly affect the choice of the 
optimal a and 0. The major increases in the convergence rate are for the Euler equations. For the 
Navier-Stokes equations it is necessary to reformulate the preconditioning matrix to account for 
the viscous effecrs. 
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3. Mesh for turbulent flow around a VKI cascade. 
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